suppressWarnings(library(gplots))
suppressWarnings(library(multcomp))
suppressWarnings(library(R2HTML))
suppressWarnings(library(lattice))

# retrieve args
Args <- commandArgs(TRUE)

#Read in data
statdata <- read.csv(Args[3], header=TRUE, sep=",")

#Copy Args
response <- Args[4]
responseTransform <- Args[5]
treatFactor <- Args[6]
equalCase <-Args[7]
unequalCase <- Args[8]
showPRPlot <- Args[9]
showNormPlot <- Args[10]
sig <- 1 - as.numeric(Args[11])

#Setup the html file and associated css file
htmlFile <- sub(".csv", ".html", Args[3]); #determine the file name of the html file
.HTML.file = htmlFile
cssFile <- "r2html.css"
cssFile <- paste("'",cssFile,"'", sep="") #need to enclose in quotes when path has spaces in it
HTMLCSS(CSSfile = cssFile)

#Output HTML header
HTML.title("<bf>SilveR Unpaired t-test Analysis", HR=1, align="left")


#Removing illegal charaters

YAxisTitle <-response
XAxisTitle <-treatFactor

for (i in 1:10)
{

# Additional characters included Aug 2010 (STB)
YAxisTitle<-sub("ivs_tilde_ivs"	,"~", YAxisTitle) 
YAxisTitle<-sub("ivs_star_ivs"	,"*", YAxisTitle) 
YAxisTitle<-sub("ivs_plus_ivs"	,"+", YAxisTitle) 

YAxisTitle<-sub("ivs_sp_ivs"	," ", YAxisTitle) 
YAxisTitle<-sub("ivs_ob_ivs"	,"(", YAxisTitle) 
YAxisTitle<-sub("ivs_cb_ivs"	,")", YAxisTitle) 
YAxisTitle<-sub("ivs_div_ivs"	,"/", YAxisTitle) 
YAxisTitle<-sub("ivs_pc_ivs"	,"%", YAxisTitle) 
YAxisTitle<-sub("ivs_hash_ivs"	,"#", YAxisTitle) 
YAxisTitle<-sub("ivs_pt_ivs"	,".", YAxisTitle) 
YAxisTitle<-sub("ivs_hyphen_ivs","-", YAxisTitle) 
YAxisTitle<-sub("ivs_at_ivs"	,"@", YAxisTitle) 
YAxisTitle<-sub("ivs_colon_ivs"	,":", YAxisTitle) 
YAxisTitle<-sub("ivs_exclam_ivs","!", YAxisTitle) 
YAxisTitle<-sub("ivs_quote_ivs"	,"`", YAxisTitle) 
YAxisTitle<-sub("ivs_pound_ivs"	,"£", YAxisTitle) 
YAxisTitle<-sub("ivs_dollar_ivs","$", YAxisTitle) 
YAxisTitle<-sub("ivs_hat_ivs"	,"^", YAxisTitle) 
YAxisTitle<-sub("ivs_amper_ivs"	,"&", YAxisTitle) 
YAxisTitle<-sub("ivs_obrace_ivs","{", YAxisTitle) 
YAxisTitle<-sub("ivs_cbrace_ivs","}", YAxisTitle) 
YAxisTitle<-sub("ivs_semi_ivs"	,";", YAxisTitle) 
YAxisTitle<-sub("ivs_pipe_ivs"	,"|", YAxisTitle) 
YAxisTitle<-sub("ivs_slash_ivs"	,"\\", YAxisTitle) 
YAxisTitle<-sub("ivs_osb_ivs"	,"[", YAxisTitle) 
YAxisTitle<-sub("ivs_csb_ivs"	,"]", YAxisTitle) 
YAxisTitle<-sub("ivs_eq_ivs"	,"=", YAxisTitle) 
YAxisTitle<-sub("ivs_lt_ivs"	,"<", YAxisTitle) 
YAxisTitle<-sub("ivs_gt_ivs"	,">", YAxisTitle) 
YAxisTitle<-sub("ivs_dblquote_ivs"	,"\"", YAxisTitle) 

# Additional characters included Aug 2010 (STB)
XAxisTitle<-sub("ivs_tilde_ivs"	,"~", XAxisTitle) 
XAxisTitle<-sub("ivs_star_ivs"	,"*", XAxisTitle) 
XAxisTitle<-sub("ivs_plus_ivs"	,"+", XAxisTitle) 

XAxisTitle<-sub("ivs_sp_ivs"	," ", XAxisTitle) 
XAxisTitle<-sub("ivs_ob_ivs"	,"(", XAxisTitle) 
XAxisTitle<-sub("ivs_cb_ivs"	,")", XAxisTitle) 
XAxisTitle<-sub("ivs_div_ivs"	,"/", XAxisTitle) 
XAxisTitle<-sub("ivs_pc_ivs"	,"%", XAxisTitle) 
XAxisTitle<-sub("ivs_hash_ivs"	,"#", XAxisTitle) 
XAxisTitle<-sub("ivs_pt_ivs"	,".", XAxisTitle) 
XAxisTitle<-sub("ivs_hyphen_ivs","-", XAxisTitle) 
XAxisTitle<-sub("ivs_at_ivs"	,"@", XAxisTitle) 
XAxisTitle<-sub("ivs_colon_ivs"	,":", XAxisTitle) 
XAxisTitle<-sub("ivs_exclam_ivs","!", XAxisTitle) 
XAxisTitle<-sub("ivs_quote_ivs"	,"`", XAxisTitle) 
XAxisTitle<-sub("ivs_pound_ivs"	,"£", XAxisTitle) 
XAxisTitle<-sub("ivs_dollar_ivs","$", XAxisTitle) 
XAxisTitle<-sub("ivs_hat_ivs"	,"^", XAxisTitle) 
XAxisTitle<-sub("ivs_amper_ivs"	,"&", XAxisTitle) 
XAxisTitle<-sub("ivs_obrace_ivs","{", XAxisTitle) 
XAxisTitle<-sub("ivs_cbrace_ivs","}", XAxisTitle) 
XAxisTitle<-sub("ivs_semi_ivs"	,";", XAxisTitle) 
XAxisTitle<-sub("ivs_pipe_ivs"	,"|", XAxisTitle) 
XAxisTitle<-sub("ivs_slash_ivs"	,"\\", XAxisTitle) 
XAxisTitle<-sub("ivs_osb_ivs"	,"[", XAxisTitle) 
XAxisTitle<-sub("ivs_csb_ivs"	,"]", XAxisTitle) 
XAxisTitle<-sub("ivs_eq_ivs"	,"=", XAxisTitle) 
XAxisTitle<-sub("ivs_lt_ivs"	,"<", XAxisTitle) 
XAxisTitle<-sub("ivs_gt_ivs"	,">", XAxisTitle)  
XAxisTitle<-sub("ivs_dblquote_ivs"	,"\"", XAxisTitle) 
}


#Generate mainEffect factor
statdata$mainEffect<-as.factor(eval(parse(text = paste("statdata$", treatFactor))))


#Warning
title<-c("Warning")
HTML.title(title, HR=2, align="left")

HTML.title("</bf> ", HR=2, align="left")
HTML.title("Warning: This module is currently under construction, care should be taken when considering the results. The results have not been verified.", HR=0, align="left")


#Response

title<-c("Response")

HTML.title(title, HR=2, align="left")
add<-paste(c("The  "), response, sep="")
add<-paste(add, " response is currently being analysed by the Unpaired t-test Analysis module.", sep="")

HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")

if (responseTransform != "None")
{
	HTML.title("<bf> ", HR=2, align="left")
}

if (responseTransform != "None")
{
	add2<-paste(c("The response has been "), responseTransform, sep="")
	add2<-paste(add2, " transformed prior to analysis.", sep="")
	HTML.title(add2, HR=0, align="left")
}


#Scatterplot
HTMLbr()
title<-c("Scatterplot of the raw data")
if(responseTransform != "None")
{
	title<-paste(title, " (on the ", sep="")
	title<-paste(title, responseTransform, sep="")
	title<-paste(title, " scale)", sep="")
}

HTML.title(title, HR=2, align="left")
scatterPlot <- sub(".html", "scatterPlot.jpg", htmlFile)
jpeg(scatterPlot)
xyplot(eval(parse(text = paste("statdata$", response)))~eval(parse(text = paste("statdata$", treatFactor))), statdata, col="black", ylab=YAxisTitle, xlab=NULL, scales=list(x=list(rot=90)) )
void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", scatterPlot), Align="centre")

HTML.title("</bf> ", HR=2, align="left")
HTML.title("</bf> Tip: Use this plot to identify possible outliers.", HR=0, align="left")


#Equal variance case

if (equalCase == "Y")
{
	HTMLbr()
	HTML.title("<bf>t-test assuming equal variances", HR=2, align="left")

	eqtest<-t.test(formula = eval(parse(text = paste("statdata$", response)))~eval(parse(text = paste("statdata$", treatFactor))), paired = FALSE, var.equal= TRUE, conf.level= sig)

	col1<-format(round(as.numeric(eqtest[1]), 3), nsmall=3, scientific=FALSE)
	col2<-format(round(as.numeric(eqtest[2]), 0), nsmall=0, scientific=FALSE)
	col3<-format(round(as.numeric(eqtest[3]), 4), nsmall=4, scientific=FALSE)

	blanq<-c(" ")
	for (i in 1:1)
	{ 
		blanq[i]=" "
	}

	eqtable<-cbind(col1,blanq,col2,blanq,col3)

	head<-c("t-statistic", " ","Degrees of Freedom"," ","p-value")
	colnames(eqtable)<-head
	rownames(eqtable)<-c("Equal variance t-test")

	for (i in 1:(dim(eqtable)[1])) 
	{
		if (eqtable[i,5]<0.001) 
		{
			eqtable[i,5]<-0.001
			eqtable[i,5]<- paste("<",eqtable[i,5])
		}
	}

	HTML(eqtable, classfirstline="second", align="left")

	add<-paste(c("Conclusion"))
	if (eqtable[i,5]<= (1-sig))
	{
		add<-paste(add, ": There is a statistically significant difference between the levels of ", sep="")
		add<-paste(add, response, " at the ", 100*(1-sig) , "% level" , sep="")
		add<-paste(add, ".", sep="")
	} else {
		add<-paste(add, ": The difference between the levels of ", sep="")
		add<-paste(add, response, sep="")
		add<-paste(add, " is not statistically significant.", sep="")
		}
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(add, HR=0, align="left")

# Means and Planned comparisons on the main effects

#Code to calculate y-axis offset (lens) in LSMeans plot
	names<-levels(eval(parse(text = paste("statdata$", treatFactor))))
	index<-1
	for (i in 1:length(names))
	{
		temp<-names[i]
		temp<-as.character(unlist(strsplit(as.character(names[i]),"")))
		lens<-length(temp)
		if (lens>index)
		{
			index <-lens
		}
	}

	mult<-glht(lm(eval(parse(text = paste("statdata$", response)))~ mainEffect, data=statdata, na.action = na.omit), linfct=mcp(mainEffect="Means"))
	meanPlot <- sub(".html", "meanplot.jpg", htmlFile)
	jpeg(meanPlot)

	tab1<-c()
	tab2<-c()
	tab3<-c()

	tabs<-confint(mult,level = sig, calpha = univariate_calpha() )
	test<-tabs$confint
	lengths<-length(unique(rownames(test)))

	for (i in 1:lengths)
	{
		tab1[i]<-test[i,1]
		tab2[i]<-test[i,2]
		tab3[i]<-test[i,3]
	}

	tab1<-format(round(tab1, 2), nsmall=2, scientific=FALSE)
	tab2<-format(round(tab2, 2), nsmall=2, scientific=FALSE)
	tab3<-format(round(tab3, 2), nsmall=2, scientific=FALSE)

	table<-c(1:lengths)
	for (i in 1:lengths)
	{
		table[i]=" "
	}

	tab <- cbind(table, tab1, table, tab2, table, tab3)

	colnames(tab)<-c(" ", "Mean", " ", paste("Lower ",(sig*100),"% CI",sep=""), " ", paste("Upper ",(sig*100),"% CI",sep=""))
	header<-c(" ", " "," "," ", " "," ")
	tables<-rbind(header, tab)
	rownames(tables)<-c("Treatment level", rownames(test))

	HTMLbr()
	CITitle2<-paste("<bf>Table of the predicted means with ",(sig*100),"% confidence intervals (assuming equal variances)",sep="")
	HTML.title(CITitle2, HR=2, align="left")

	HTML(tables, classfirstline="second", align="left")

#Testing the degrees of freedom

	if (eqtest[2]<5)
	{
		HTMLbr()
		HTML.title("<bf>Warning", HR=2, align="left")
	
		add<-c("Unfortunately the residual degrees of freedom are low (less than 5). This may make the estimation of the underlying variability, and hence the results of the statistical tests, unreliable.")
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(add, HR=0, align="left")
		print(add)
	}
}

#Unequal variance case

if (unequalCase == "Y")
{
	HTMLbr()
	HTML.title("<bf>t-test assuming unequal variances", HR=2, align="left")

	add<-c("The analysis presented is Welch's t-test for unequal variances.")
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(add, HR=0, align="left")

	eqtest<-t.test(formula = eval(parse(text = paste("statdata$", response)))~eval(parse(text = paste("statdata$", treatFactor))), paired = FALSE, var.equal= FALSE, conf.level= sig)

	col1<-format(round(as.numeric(eqtest[1]), 3), nsmall=3, scientific=FALSE)
	col2<-format(round(as.numeric(eqtest[2]), 2), nsmall=2, scientific=FALSE)
	col3<-format(round(as.numeric(eqtest[3]), 4), nsmall=4, scientific=FALSE)

	blanq<-c(" ")
	for (i in 1:1)
	{ 
		blanq[i]=" "
	}

	eqtable<-cbind(col1,blanq,col2,blanq,col3)

	head<-c("t-statistic", " ","Degrees of Freedom"," ","p-value")
	colnames(eqtable)<-head
	rownames(eqtable)<-c("Unequal variance t-test")

	for (i in 1:(dim(eqtable)[1])) 
	{
		if (eqtable[i,5]<0.001) 
		{
			eqtable[i,5]<-0.001
			eqtable[i,5]<- paste("<",eqtable[i,5])
		}
	}

	HTML(eqtable, classfirstline="second", align="left")

	add<-paste(c("Conclusion"))
	if (eqtable[i,5]<= (1-sig))
	{
		add<-paste(add, ": There is a statistically significant difference between the levels of ", sep="")
		add<-paste(add, response, " at the ", 100*(1-sig) , "% level" , sep="")
		add<-paste(add, ".", sep="")
	} else {
		add<-paste(add, ": The difference between the levels of ", sep="")
		add<-paste(add, response, sep="")
		add<-paste(add, " is not statistically significant.", sep="")
		}
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(add, HR=0, align="left")

# Means and Planned comparisons on the main effects
	tablenames<-c(levels(as.factor(statdata$mainEffect)))
	length <- length(unique(levels(as.factor(statdata$mainEffect))))
	table<-c(1:length)
	for (i in 1:length)
	{
		table[i]=" "
	}
	vectorUCI <-c(1:length)
	vectorLCI <-c(1:length)
	vectormean <-c(1:length)
	for (i in 1:length)
	{
		sub<-subset(statdata, statdata$mainEffect == unique(levels(as.factor(statdata$mainEffect)))[i])
		sub2<-data.frame(sub)
		tempy<-na.omit(eval(parse(text = paste("sub2$", response))))
		vectormean[i] = mean(eval(parse(text = paste("sub2$", response))), na.rm=TRUE)
		vectorLCI[i]= mean(eval(parse(text = paste("sub2$", response))), na.rm=TRUE)-qt(1-(1-sig)/2, (length(tempy)-1))*sd(eval(parse(text = paste("sub2$", response))), na.rm=TRUE) / (length(tempy))**(0.5)
		vectorUCI[i]= mean(eval(parse(text = paste("sub2$", response))), na.rm=TRUE)+qt(1-(1-sig)/2, (length(tempy)-1))*sd(eval(parse(text = paste("sub2$", response))), na.rm=TRUE) / (length(tempy))**(0.5)
	}

	vectormean<-format(round(vectormean, 2), nsmall=2, scientific=FALSE)
	table2<-cbind(table,vectormean)
	vectorLCI<-format(round(vectorLCI, 2), nsmall=2, scientific=FALSE)
	table2<-cbind(table2,table)
	table2<-cbind(table2,vectorLCI)
	vectorUCI<-format(round(vectorUCI, 2), nsmall=2, scientific=FALSE)
	table2<-cbind(table2,table)
	table2<-cbind(table2,vectorUCI)

	header<-c(" "," "," ", " ", " ", " ")
	temp6<-c("Mean")
	temp7<-c(" ")

	CIlow<-paste("Lower ", 100*(sig), sep="")
	CIlow<-paste(CIlow , "% CI", sep="")
	hed12<-c(CIlow)

	temp6<-cbind(temp7, temp6,temp7)
	temp6<-cbind(temp6,hed12)
	temp6<-cbind(temp6,temp7)

	CIhigh<-paste("Upper ", 100*(sig), sep="")
	CIhigh<-paste(CIhigh, "% CI", sep="")
	hed13<-c(CIhigh)
	temp6<-cbind(temp6,hed13)
	table2<-rbind(header, table2)

	colnames(table2)<-temp6
	rownms<-c("Treatment level")
	for(i in 1:length)
	{
		rownms[i+1]<-levels(as.factor(statdata$mainEffect))[i]
	}
	
	row.names(table2)<-rownms

	HTMLbr()
	CITitle2<-paste("<bf>Table of the predicted means with ",(sig*100),"% confidence intervals (assuming unequal variances)",sep="")
	HTML.title(CITitle2, HR=2, align="left")

	HTML(table2, , align="left" , classfirstline="second")
	HTML.title("</bf> ", HR=2, align="left")

#Testing the degrees of freedom
	if (eqtest[2]<5)
	{
		HTMLbr()
		HTML.title("<bf>Warning", HR=2, align="left")
	
		add<-c("Unfortunately the residual degrees of freedom are low (less than 5). This may make the estimation of the underlying variability, and hence the results of the statistical tests, unreliable.")
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(add, HR=0, align="left")
		print(add)
	}
}


#Diagnostic plot titles
if(showPRPlot=="Y")
{
	HTMLbr()
	HTML.title("<bf>Diagnostic plots", HR=2, align="left")
} else {
	if(showNormPlot=="Y")
	{
		HTMLbr()
		HTML.title("<bf>Diagnostic plots", HR=2, align="left")
	}
}

#Residual plots
if(showPRPlot=="Y")
{
	threewayfull<-lm(eval(parse(text = paste("statdata$", response)))~ mainEffect, data=statdata, na.action = na.omit)

	residualPlot <- sub(".html", "residualplot.jpg", htmlFile)
	jpeg(residualPlot)

	residplot<-data.frame(cbind(predict(threewayfull),rstudent(threewayfull)))
	colnames(residplot)<-c("Predicted","Standardized_residuals")
	plot(Standardized_residuals~Predicted,data=residplot, ylab="Externally Studentized residuals", xlab="Predicted values",main="Predicted vs. residuals plot")
	abline(a=0, b=0)
	abline(a=2, b=0, col="red",lty=3) 
	abline(a=3, b=0, col="red")
	abline(a=-2, b=0, col="red",lty=3)
	abline(a=-3, b=0, col="red")

	abline(h=0, col="red", lty="dotted")
	void<-HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", residualPlot), Align="centre")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>Tip: On this plot look to see if the spread of the points increases as the predicted values increase. If so the response may need transforming or the unequal variance assumption selected.", HR=0, align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>Tip: Any observation with a residual less than -3 or greater than 3 (SD) should be investigated as a possible outlier.", HR=0, align="left")
}

#Normality plots
if(showNormPlot=="Y")
{
	threewayfull<-lm(eval(parse(text = paste("statdata$", response)))~ mainEffect, data=statdata, na.action = na.omit)

	HTMLbr()
	normPlot <- sub(".html", "normplot.jpg", htmlFile)
	jpeg(normPlot)
	qqnorm(resid(threewayfull) ,main="Normal probability plot")
	qqline(resid(threewayfull), col="red", lty="dotted")
	void<-HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", normPlot), Align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>Tip: Check that the points lie along the dotted line. If not then the data may be non-normally distributed.", HR=0, align="left")
}

#References

if(unequalCase == "Y")
{
HTMLbr()
HTML.title("<bf>Statistical references", HR=2, align="left")

HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf> Welch BL. (1947). The generalization of Student's problem when several different population variances are involved. Biometrika, 34(1-2), 28-35.", HR=0, align="left")
}

HTMLbr()
HTML.title("<bf>R references", HR=2, align="left")

HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>   R Development Core Team (2008). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.", HR=0, align="left")

#lattice
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf> 
	Deepayan Sarkar (2009). lattice: Lattice Graphics. R package version 0.17-22. http://CRAN.R-project.org/package=lattice
	", HR=0, align="left")

#gplots
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
 	Gregory R. Warnes. Includes R source code and/or documentation contributed by (in alphabetical order): Ben Bolker, Lodewijk Bonebakker, Robert Gentleman, Wolfgang Huber, Andy Liaw, Thomas Lumley, Martin Maechler, Arni Magnusson, Steffen Moeller, Marc Schwartz and Bill Venables (2009). gplots: Various R programming tools for plotting data. R package version 2.7.1. http://CRAN.R-project.org/package=gplots
	", HR=0, align="left")

#R2HTML
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
	Lecoutre, Eric (2003). The R2HTML Package. R News, Vol 3. N. 3, Vienna, Austria.
	", HR=0, align="left")

#multcomp
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>    
	Torsten Hothorn, Frank Bretz and Peter Westfall with contributions by Richard M. Heiberger (2007). multcomp: Simultaneous Inference for	General Linear Hypotheses. R package version 0.991-8.
	", HR=0, align="left")